Implementation of higher-order absorbing boundary 
conditions for the Einstein equations 



Oliver Rinne^^ Luisa T. Buchman'^^^, Mark A. ScheeF 
and Harald P. PfeifFer^ 

^Department of Applied Mathematics and Theoretical Physics, Centre for 
Mathematical Sciences, Wilberforce Road, Cambridge CBS OWA, UK 
^King's College, Cambridge CB2 1ST, UK 

■^Theoretical Astrophysics 130-33, California Institute of Technology, Pasadena, 
CA 91125, USA 

''Center for Relativity, University of Texas at Austin, Austin, TX 78712, USA 



Abstract. We present an implementation of absorbing boundary conditions 
for the Einstein equations based on the recent work of Buchman and Sarbach. 
In this paper, we assume that spacetime may be linearized about Minkowski 
space close to the outer boundary, which is taken to be a coordinate sphere. We 
reformulate the boundary conditions as conditions on the gauge-invariant Regge- 
Wheeler-Zerilli scalars. Higher-order radial derivatives are eliminated by rewriting 
the boundary conditions as a system of ODEs for a set of auxiliary variables 
intrinsic to the boundary. From these we construct boundary data for a set of 
well-posed constraint-preserving boundary conditions for the Einstein equations 
in a first-order generalized harmonic formulation. This construction has direct 
applications to outer boundary conditions in simulations of isolated systems (e.g., 
binary black holes) as well as to the problem of Cauchy-perturbative matching. 
As a test problem for our numerical implementation, we consider linearized 
multipolar gravitational waves in TT gauge, with angular momentum numbers 
1 = 1 (Teukolsky waves), 3 and 4. We demonstrate that the perfectly absorbing 
boundary condition Bi^ of order L = £ yields no spurious reflections to linear order 
in perturbation theory. This is in contrast to the lower-order absorbing boundary 
conditions Bl with L < £, which include the widely used freezing- 'to boundary 
condition that imposes the vanishing of the Newman-Penrose scalar 'I'o. 



PACS numbers: 04.25. D-, 02.60.Lj, 04.25.-g 

1. Introduction 



Many situations of astropliysical interest can be described to good approximation as an 
isolated system: an asymptotically flat spacetime containing a compact self-gravitating 
source. The study of such systems requires the solution of Einstein's equations on 
an unbounded domain. One of the major problems in numerical relativity is how 
to accomplish this with finite computer resources. The most common approach is 
based on the Cauchy formulation of general relativity. Here the spatial computational 
domain is truncated at a finite distance and boundary conditions (BCs) are imposed 
at this artificial timelike boundary. These BCs must satisfy a number of requirements, 
the most important ones being that (i) the BCs must be compatible with the constraint 
equations that hold within the spatial slices, i.e. the BCs must be constraint preserving, 
(ii) the resulting initial-boundary value problem (IBVP) must be well posed, and 
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(iii) the BCs should minimize spurious reflections of gravitational radiation, i.e. the 
BCs should be absorbing. (We call BCs that completely eliminate such reflections 
perfectly absorbing.) It is important here to distinguish between spurious reflections 
on the one hand, and non-spurious reflections such as the backscatter off a curved 
background spacetime on the other hand. Perfectly absorbing BCs are defined to 
be exactly satisfied by the general retarded solution at the outer boundary, i.e. they 
only eliminate the spurious reflections but preserve non-spurious reflections such as 
backscatter. While considerable progress has been made on the first two requirements 
above, i.e. constraint preservation and well posedness, the third one on the absorbing 
properties of the BCs has not been addressed to the same extent. It is however of 
prime importance if accurate approximations to the gravitational radiation emitted 
by the source are to be computed, as required for instance by gravitational wave data 
analysis. 

Theoretical progress on the construction of absorbing BCs for the Einstein field 
equations was made by Buchman and Sarbach in [HIl]. A hierarchy of local BCs Bl 
was proposed that is perfectly absorbing for all radiation with angular momentum 
numbers £ ^ L, where L is an arbitrary given number. These BCs were obtained 
by studying solutions to the Bianchi equations describing linearized gravitational 
waves and finding a condition on the outgoing solutions that was satisfied exactly. 
It turns out that these BCs are a generalization to the Einstein equations of the well- 
known Bayliss-Turkel BCs [3] for the scalar wave equation. Initially [T], the BCs were 
formulated for a flat background spacetime; first-order corrections dealing with the 
spacetime curvature and the backscatter on a Schwarzschild black hole background 
were included in [2]. 

The objective of the present paper is to reformulate the Buchman-Sarbach BCs Bl 
in such a way that they can be incorporated into a full set of constraint-preserving BCs 
for the Einstein equations, and to implement and test such BCs. Before describing our 
approach in more detail, we mention a number of previous works on and alternatives 
to absorbing BCs. 

In [3], the Bayliss-Turkel BCs were implemented up to quadrupolar order (£ = 2) 
for the scalar wave equation. The reformulation of the BCs and the numerical method 
used in that paper are similar to ours, but we note that we treat gravitational rather 
than scalar waves. A different approach to absorbing BCs is based on fast-converging 
series expansions of exact nonlocal BCs [5]. This was applied to the construction of 
exact absorbing BCs for the Regge- Wheeler and Zcrilli (RWZ) equations in a series of 
papers by Lau [HllTllHj- The RWZ equations describe linear gravitational perturbations 
about a Schwarzschild black hole, and they play a central role in our approach as 
well. A general framework for the construction of absorbing boundaries is Cauchy- 
perturbative matching (CPM) [9l [TOl HU [12]. Here one matches solutions of the 
Einstein equations in the interior of a compact domain to solutions of the linearized 
equations in the exterior, represented by, for example, the RWZ equations. As we shall 
point out below, this approach is again closely related to ours. In addition, instead of 
matching to solutions of linearized gravity, one can match to an outer module solving 
the Einstein equations in the characteristic formulation, which is particularly well 
suited to the extraction of gravitational radiation (see [13] for a review article) . Finally, 
a very promising method consists in solving the Einstein equations on hyperboloidal 
slices that are everywhere spacelike but that asymptotically approach null infinity (see 
[14] for a review article). 

Our approach is based on the Buchman-Sarbach BCs, and in the following we 
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describe the main idea of our algorithm and the organization of this paper in more 
detail. We assume that close to the outer boundary, the spacetime metric can be 
described by linear perturbations of Minkowski spacetime in the standard coordinates. 
This is usually a good approximation if the outer boundary is placed sufficiently far out 
(we intend to generalize our work to a Schwarzschild background in the future). Note 
that on a flat background (unlike on a curved background where there is backscatter), 
the general retarded solution is purely outgoing, so in this case perfectly absorbing BCs 
should eliminate precisely the ingoing solution. We furthermore assume that the outer 
boundary is a sphere of constant coordinate radius. The original Buchman-Sarbach 
conditions involve higher-order radial derivatives of the Newman-Penrose scalar 5*0, 
making it a non-trivial task to implement them numerically. Our strategy is to work 
with the RWZ scalars instead of ^'o because on a flat-space background, the Buchman- 
Sarbach BCs on the RWZ scalars are precisely the same as the Bayliss-Turkel BCs for 
the scalar wave equation (section lTTI) . Additionally, since the RWZ scalars obey closed 
wavelike evolution equations (the RWZ equations), we can draw on previous work in 
computational mathematics that successfully implements higher-order absorbing BCs 
for the scalar wave equation. Following [151 HSl HIIj we eliminate the higher-order 
radial derivatives by rewriting the BCs as a system of ordinary differential equations 
(ODEs) for a set of auxiliary variables that need only be deflned at the boundary 
(section l2.2p . This auxiliary system is completely equivalent to the absorbing BCs 
imposed on the RWZ scalars. 

In the interior of the computational domain, the full (nonlinear) Einstein 
equations are evolved. Hence we must transfer information between them and the 
auxiliary system at the boundary. The RWZ scalars are computed from the spacetime 
metric in section 12.31 following the gauge-invariant treatment of Sarbach and Tiglio 
|18| . In turn, we must provide boundary data for the incoming characteristic fields 
of the Einstein equations (section 12. 4p . We use a flrst-order formulation of these 
equations in generalized harmonic gauge (see [Hj and references therein). In |20j . 
a set of constraint-preserving and well-posed BCs in Sommerfeld form was proposed 
for the harmonic Einstein equations. These BCs contain certain free boundary data, 
and only very special values of these data will yield BCs that are also absorbing. 
The novelty of our approach lies in the construction of such absorbing boundary data 
from the auxiliary system at the boundary (section 12. 5|) . We stress here that these 
data could equally well be obtained from exterior solutions of the RWZ equations in 
the CPM approach, so that our work can also be viewed as an explicit prescription 
for constructing BCs for the generalized harmonic Einstein equations from such 
perturbative solutions. 

In summary, our algorithm consists in three basic steps: (i) extraction of the 
RWZ scalars from the spacetime metric at the boundary, (ii) evolution of the auxiliary 
variables at the boundary, and (iii) construction of absorbing boundary data for the 
Einstein equations from the auxiliary variables. 

We evolve the generalized harmonic Einstein equations using a pseudospectral 
collocation method, the Caltech-Cornell Spectral Einstein Code (SpEC). Some details 
of the numerical method relevant to the present work are described in section [3.31 

As a first test problem of our implementation, we consider exact solutions of 
linearized gravity, multipolar gravitational waves in transverse-traceless (TT) gauge 
(section 13. ip . The quadrupolar (£ = 2) solution was first given in explicit form by 
Teukolsky [21] , and it has recently been generalized to arbitrary £ in [22| . The BC BL=e 
is perfectly absorbing for these solutions. We set initial data for £ = 2,3,4 (section 
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13. 2p and evolve them with this BC, B^^i, comparing the numerically extracted RWZ 
scalars with their analytical counterparts (section 13. 5p . We show that the difference 
between the two decays (at least) quadratically with amplitude so that there are no 
spurious reflections to linear order in perturbation theory, as expected. In contrast, 
BL<e are shown to cause reflections at leading (linear) order. The L = 1 case, Si, is 
equivalent to imposing the vanishing of the Newman-Penrose scalar ^'o at the outer 
boundary, a BC that is often used in numerical relativity [23l [24l [25 l [27l [28 l [29] 
and that is referred to as the freezing-5'o BC. There is a residual gauge freedom 
at the boundary in generalized harmonic gauge (section 13. 4p . and we show that our 
numerical results are insensitive to the particular choice of gauge boundary data. We 
also compute approximate reflection coefflcients from our numerical evolutions and 
compare them with the theoretical predictions of [T] (section [3.6p . 
We conclude and give an outlook onto future work in section IH 

2. Formulation of the boundary conditions 

In this section, we derive our reformulation of the higher-order absorbing BCs Bl 
proposed in [U [2]. We begin by showing how they can be expressed as BCs on 
the RWZ scalars (section 12. 1|) . Radial derivatives are eliminated by introducing a 
system of auxiliary ODEs at the boundary (section l2.2p . Next we explain how the 
RWZ scalars are extracted from the perturbed spacetime metric, following the gauge- 
invariant treatment of Sarbach and Tiglio [18] (section 12. 3[) . Finally we describe the 
BCs that we impose on the Einstein equations in a first-order generalized harmonic 
formulation of the Einstein equations (section 12. 4p . Boundary data are constructed 
from the auxiliary system that correspond precisely to the desired absorbing BC 
(section 12. Sp . 

2.1. Absorbing boundary conditions for the RWZ scalars 

Although our approach is not limited to this case, we assume in this paper that 
spacetime near the outer boundary can be described by linear perturbations of a flat 
background spacetime in standard Minkowski coordinates. The spacetime metric gap 
is written as 

gap ^ gap + 5gaf3. (1) 

We assume that the background metric gap is a direct product 

g = gabdx''dx'' + r^gAsdx'^dx^ , (2) 

where g = —dt^ + dr^ is the standard Minkowski metric on a 2-manifold M and 
g = d9^ + sin^ 9 d(jp' is the standard metric on the 2-sphere. Throughout this paper, 
Greek indices a, /?, . . . are spacetime indices, lower-case Latin indices a, 6, . . . range 
over t and r, and upper-case Latin indices A, B, . . . range over 9 and (f>. The covariant 
derivative compatible with the metric g (g) will be denoted by V (V) and the volume 
element by e^b {^ab)- The covariant derivative V is sometimes also denoted by a 
vertical bar (|). 

The hierarchy of absorbing BCs Bl proposed in [1] are, for a flat-space 
background, 

Bl: [r'{dt + dr)]'^-\r'^o) = 0, (3) 
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where = denotes equahty at the boundary. Here the Newman-Penrose scalar 
is evaluated for a null tetrad (P, A:", to", m") adapted to the background spacetime. 
The BCs ^ are perfectly absorbing for all perturbations with angular momentum 
numbers £ ^ L. We note that for L ~ 1, ^ reduces to the often-used freezing- 
condition, 9t*o = 0[2l[24l[25l[271[26l[28l[29]. (We have not included the outer time 
derivative in Q ; the purpose of that time derivative is to address a static background 
contribution to ^'o, but this contribution vanishes in flat space.) 

Instead of ^^o, we choose to work with the RWZ scalars describing gauge- 
invariant gravitational perturbations about Schwarzschild spacetime (see |18] and 
references therein; of course, this includes our assumption of a flat-space background 
as a special case). The superscript (±) refers to the parity of the perturbations: 
(-I-) for even and (— ) for odd parity. The indices £m refer to a spherical harmonic 
decomposition and will usually be suppressed in the following. We use the RWZ 
scalars because they have the advantage that they obey a closed evolution equation, 
the RWZ equation, which in flat space reads 



= 0. (4) 



This equation arises from the scalar wave equation on /r after a decomposition 
into spherical harmonics; it is also known as the Euler-Poisson-Darboux equation |30j . 

The relation between and the perturbations J^'o of is provided by 

equation 22 in [2], which for a flat-space background reduces to 

5*0 = rZ^VaVb[r($(+) + i$(-))]TO'^TO'°VcVi)r. (5) 

Here Y are the standard scalar spherical harmonics. We have suppressed indices £m 
on Y and 'I'^^-', which are being summed over. 

Equation ^ allows us to translate the BCs ^ on *o into BCs on 

P(9t+a,)]^+i$ = 0, (6) 

which holds for both parities (and hence we suppress the superscript (±)). We note 
that ([6|) are the well-known Bayliss-Turkel conditions [3] for the scalar wave equation 
on $/?■, equation (j4]). 

2.2. The auxiliary system at the boundary 

The BCs ([6]), which arc equivalent to the Bl conditions ([3]), are difficult to implement 
numerically because they contain higher-order radial derivatives. We follow a 
procedure developed for the scalar wave equation in [151 1161 117] in order to eliminate 
these derivatives. A set of auxiliary variables is introduced, 

Wk^r-'^^'^+'^lr^dt+dr)]"^, (7) 

where again the parity (±) and the indices £m are suppressed. Hence these auxiliary 
variables obey the recursion relation 

2fc + 1 \ 

dt+dr^ Wk = wk+i- (8) 



r 

Using the wave equation (|4]) and induction over fc, we can prove the identity [16| 

dt~dr^ -) Wk = ^[-£{£+1) + k{k - l)]wk-i. (9) 
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Equations ^ and Q can be combined to eliminate the radial derivatives, 

{^t + ^) w'fc = ^i-^i^ + 1) + Hk - l)]wk-i + ^Wk+i. (10) 
The BC (III) is equivalent to 

WL+1=0, (11) 



which closes the system of ODEs ([T0|) . We integrate (|T0|) on the boundary for 
1 < fc < L, substituting pT|) and wq — $/r. 

^.5. Extraction of the RWZ scalars 

The RWZ scalars need to be computed from the spacetime metric. We follow the 
gauge-invariant treatment of [18] . restricted to a Minkowski background in standard 
coordinates. The starting point is a decomposition of the metric perturbations with 
respect to scalar, vector and tensor spherical harmonics. The basis harmonics are 
defined by 

Ya = VaY, 5^ = e^A>B, 

Yab = [V(ArB)]T^ - ^iA^B}Y + + lygABY, Sab = ^i^aSb), (12) 
where TF denotes the tracefree part. The two parities are treated separately. 

2.3.1. Odd parity. Odd-parity perturbations of the spacetime metric are decomposed 
as 

^gAb = hbSA, (13) 
6gAB=2kSAB- (14) 
From the amplitudes ha and k, we construct the gauge- invariant potential 

h^r^-ha-r^Val^^y (15) 

i.e. 

hi"'^^ ^ht-k (16) 

(17) 

Here and in the following, a dot (prime) denotes partial differentiation with respect 
to t (r). The Regge- Wheeler scalar is defined as 



(18) 

where A = {£ — l){£ + 2). Equation psp is valid for i ^ 2; the special case 
i = 1 corresponds to a non-dynamical degree of freedom (variation of the background 
angular momentum) that is not needed in our treatment. 
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2.3.2. Even parity. Even-parity perturbations are decomposed as 

Sgab =HabY, (19) 
SgAb ^QbYA, (20) 
SgAB = r^KsABY + GYab)- (21) 

We define the gauge parameters 

Pa=Qa- yCla, (22) 

i.e. 

Pt^Qt- (23) 

Pr = Qr - (24) 

The Zcrilli 1-form is given by 

Za = Hatr^' - rA'ia - ye + l)rG\a + rl'c^ab + 2rVb\aP'' (25) 
with ujab = 2p[6|a] and Va = r^"-/r, i.e. 

Zt ^ Htr-rk ~ y£+l)rG+Pr~Pt, (26) 

Zr = Hrr - rK' ~ + l)rG" - -p^. (27) 

r 



The gauge-invariant potential is 

^(inv) ^ ^ ^ _ ^^Ihp^ ^ ^ ^ ^ _ ^p^^ (28) 

Finally, we obtain the Zerilli scalar 

= -^^(2rl"Z„ + A/.(-)) ^ -^^(2^. + A/.(-)). (29) 

Again, this formula is valid for ^ 2; the non-dynamical cases ^ = (variation of the 
background mass) and I ~ \ (pure gauge) are not needed. 

2.^. Boundary conditions for the generalized harmonic Einstein equations 

Next we describe the BCs that we impose on the actual Einstein equations that are 
being solved in the interior of the domain. We consider a first-order formulation of 
these equations in generalized harmonic coordinates (see [19| and references therein). 
Generalized harmonic spacetime coordinates x" are defined by 

□x" = -/T";37 = H", (30) 

where □ is the scalar d'Alembert operator of the spacetime metric gap, ^"p-y are its 
Christoffel symbols, and H" are freely specifiable gauge source functions. The evolved 
variables are g^p and their first derivatives 

IIq^ = - f^d^gap, (31) 

^iap = digap, (32) 

where t" is the future-directed unit timelike normal to the t — const slices. Greek 
indices a, (3, . . . are spacetime indices and Latin indices i, j, . . . from the middle of the 
alphabet are spatial. Boundary data must be specified for the incoming fields at the 
boundary, 

= n^/J - n^^taP - l2gap, (33) 
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where is the outward-pointing unit spatial normal to the boundary on the t = const 
slices, t^riQ, = 0. The parameter 72 appears due to the addition of constraint-damping 
terms to the evolution equations; let us also define 

"i/3 = + 725a/3 = ^af3 " n'^^iap. (34) 

Our BCs arc of a similar type as those considered in [20j but with boundary data 
derived from our auxiliary system at the boundary as described in the following. 

The incoming fields are split into three different pieces by three projection 
operators P*"^, and (referring to constraint, physical and gauge BCs). In order 
to define them, we introduce the null vectors = {t" +n")/y/2 a,nd k" = (P— n")/\/2 
and the spatial metric induced on the boundary. Pap = gap + tatp — UaUp. The 
projection operators are given by 

pc-ys ^ 2k^^6a^^ - kag''\ (35) 

Pip''' = Pc'^Pp' - hPc'pP''', (36) 

PG7« =Sa^l^. (37) 

The complements of the kernels of these three projection operators are disjoint and 
together they span the space of all symmetric 2-tensors. 

We impose constraint-preserving BCs by requiring the generalized harmonic gauge 
constraint ([30)1 to hold at the boundary. This can be written in the form 

pI'"K's=pS, (38) 

where Fa is a function of outgoing and zero-speed characteristic fields and gauge 
source functions, sec equation 32 of [27]. The rank-2 projection operator P^ describes 
the BCs on the two gravitational degrees of freedom. They take the form 

= Fip- (39) 

Here F^^ arc boundary data that will be constructed in the following subsection 
from our auxiliary variables at the boundary; thus these "physical" boundary data 
implement the absorbing BC Bl that we want to impose. The remaining BCs 

pG7^^;--FG (40) 

are related to the residual gauge freedom within the generalized harmonic gauge: 
we can still perform infinitesimal coordinate transformations of the coordinates 
— ^ x" + ^" provided that obeys the wave equation. The F^ are free data 
that will be specified so that our BCs are compatible with the exact solution that 
we want to reproduce. In realistic simulations without an exact solution at hand, we 
would take the F^ to vanish. 

The BCs dSHll-dlOl) were proven in [20l EI] to yield a well-posed IBVP for the 
Einstein equations in harmonic gauge. This result cannot be directly applied to our 
formulation because (i) whereas the second-order Einstein equations are considered in 
[20l[3T], our evolution system [19] is a first-order reduction thereof, and (ii) rather than 
being a priori given functions of time, the boundary data F^p in (|39|1 depend implicitly 
on the dynamical fields, as we shall see in the following subsection. Nevertheless, it 
has been shown in [32 that the absorbing BCs Bl are well posed for the second-order 
Einstein equations in harmonic gauge at least in the high-frequency limit. Therefore, 
it seems likely that our method leads to a well-posed IBVP; a proof is beyond the 
scope of this paper. 
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2.5. Construction of absorbing boundary data 

Finally, we show how to construct boundary data F^^^ in ((551) that correspond to B^. 
Recall that we are assuming that the fields can be linearized about flat space close to 
the outer boundary, which is taken to be a sphere r = const. Hence we have t" = (5t" 
and n" = 5^" for the background, and the incoming fields (without the 72 term) in 
([34]) read 

u];^3 = -^9t+dr)gap. (41) 
The only non-vanishing components of the "physical" BCs (|39p are the angular 
components, a[3 = AB. We find 

FIb = Pab'^K's - -^'(^* + dr){r-HgYB), (42) 
where TF denotes the trace-free part with respect to the metric g on the 2-sphere. 
We now derive expressions for the right-hand side of the above equation in the RWZ 
formalism that involve the auxiliary variables at the boundary. 

2.5.1. Odd parity. From p6)) and ((T7)) we obtain 

^(inv) ^ ^(inv) ^ j^^ +hr-r^ {dt + 5,) (r'^ fc) . (43) 

On the other hand [18], the gauge-invariant potential is related to the Regge- Wheeler 
scalar via /if^'"^) = ?d(r$^~)), where ? denotes the Hodge dual with respect to g. Thus, 

/if""^ + h^^^'^ = -[dt + 9.)(r$(-)) = -r^u^"^ - rw[-\ (44) 

where we have substituted the definition ([7]) of the auxiliary variables. Combining 
(03]) and dH]), 

r'^{dt + dr){r-^k) = ht + K + r"^ w[~'^ + rw[~\ (45) 
Finally (j45)) is substituted in the expression (1421) for the boundary data (noting ([M]) ). 
FIb = - r^{dt + dr){r-HgYB) = "2^^ [{dt + dr){r-^k)\ Sab 



ht + hr + r'^Wi + rwQ 



Sab- (46) 



Here we see clearly how the information from the auxiliary system that implements 
Bl enters the boundary data for the Einstein equations. We note that Sab — for 
£ ^ 1, which justifies our neglecting this special case. 



2.5.2. Even parity. From (p6)) and (|27)) wc obtain 

Zt+Zr = -^r{dr+dr)G-r{dt+dr)K + Qr-Q't-^Qr + Htr+Hrr.{^7) 

Using the definitions ^ dC = Z and = C/A, wc find the alternate expression 

Zt + Zr = \{dt + ^.O^'^' = ArwJ+^ , (48) 

where we have substituted the definition ([7]) of the auxiliary variables. Combining 
67]) and dlH]), 

r'(ft + dr)G = -2r'«;^+' - f r2(5t + dr)K + ^{Qr-Q't- iQr + Htr + Hrr). (49) 
Finally (|^9)) is substituted in the expression (|42p for the boundary data (noting (pi]) ). 
FIb = - r\dt + dr){r-HgYB) = - V^dt + dr)G] Yab 

= 2r^w[+^ - Ir^idt + dr)K + ^(Q, _ - + + H„)\ Yab- (50) 

We note that Yab = for ^ < 1, which justifies our neglecting this special case. 
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2.5.3. Application to Cauchy-perturbative matching. Wc remark that the auxihary 
variables w\^^ {k = 0, 1) appearing on the right-hand sides of (|46|) and ([50]) could 
equally well be computed from a given exterior solution of the RWZ equations, 
using the definition ([7]) of the auxiliary variables. Hence we have shown how to 
construct BCs for the generalized harmonic Einstein equations in the CPM approach. 

3. Numerical tests 

In this section, we present numerical tests of our formulation of higher-order absorbing 
BCs Bl derived in section [21 

We evolve initial data representing exact wavelike solutions of the linearized 
Einstein equations using a fully nonlinear evolution code that implements Bl- We 
extract the RWZ scalars at the boundary and compare them to the corresponding 
exact linearized solutions in order to assess the amount of spurious reflections from 
the boundary. Because we evolve the full nonlinear equations but our initial data and 
the solutions with which we compare are based on linearized theory, our comparisons 
will disagree because of nonlinear effects; we repeat each run at several different 
amplitudes in order to separate these nonlinear effects from the boundary reflections 
we are studying. 

In section 13. li we describe exact solutions of linearized gravity representing 
multipolar (with angular momentum number £) gravitational waves in TT gauge. We 
compute the RWZ scalars and demonstrate that BL=e is perfectly absorbing for these 
solutions. In section [3T2| we describe our use of these solutions to set up initial data for 
numerical evolution. After briefly explaining our pseudo-spectral numerical method 
(section [3. 3|) and our treatment of gauge BCs (section [3. 4p . we discuss our numerical 
results in section [331 Finally, in section [3?6l we estimate reflection coefficients from 
our numerical evolutions with the various BCs and compare them with the predictions 

of n. 

3.1. Multipolar gravitational wave solutions 

Tcukolsky [21] gave the explicit form of the metric for a quadrupolar {£ = 2) 
gravitational wave linearized about flat space in the TT gauge. This solution has 
recently been generalized to arbitrary angular momentum number £ in [22]. The 
solutions are characterized by freely specifiable mode functions Fout(''' — t), Fin{r + 1) 
for even parity and G'out('' — t), Gmir + t) for odd parity. They describe outgoing 
and incoming waves, respectively. There is an independent solution for each m, 
-£i^mi^£. 

We note that whereas the solution for the metric is real, the spherical harmonics 
are complex for m ^ 0. As a consequence, the amplitudes appearing in the spherical 
harmonic decomposition, equations ([T5)) - P^ and P^ - pT|) . as well as the RWZ 
scalars obey the reality conditions 

^im - i-iy'ue^-ra, (51) 

where ★ denotes complex conjugation. More precisely, let us consider the solution 
for a single mode {£,m) = {£,m). Then we may write for any of the aforementioned 
quantities 

uim = cemu{t,r), (52) 
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where u{t, r) is a real function and the Cf,„ are complex coefficients given by 

if m ^ : = '^e-,h ^ h Qm = otherwise, (53) 

if m < : = ^, C| = - 5, Qm = otherwise. (54) 

In the following, we only display the real functions u{t,r). The £ = 2 solutions in 
[2T1[22] differ by an overall constant factor; we use the normalization of [22] throughout. 

From the explicit form of the metric, equations 4 and 5 in [52], we read off the 
amplitudes 

Htt ~ Htr = 0, Hrr = A, Qt ^ 0, 

Qr = rB, K = -\A, G = C (55) 

for even parity and 

ht = 0, hr = rK, k = \r'^L (56) 

for odd parity, where the quantities B, C and K, L on the right-hand sides are the 
radial functions defined in equations 8 and 9 in [22j . 

Next we compute the RWZ scalars as described in section 12.31 For brevity, we 
focus on ^ = 2. In terms of the mode functions, the outgoing solution is found to be 



^(+)^^(4)_££^^££out^ (57) 
$(-)=Gi^J^_^ + ^, (58) 



where F('=)(x) = d^F{x)ldx^. It can be verified easily that the RWZ scalars given in 
(|57|) and ^ obey the RWZ equation (g]). 

The auxiliary variables ([T]) for the outgoing solution are found to be 

„W_3^;||_6^ .W.^, .W=0, (59) 



and similarly for odd parity. We see that the BC B2 corresponding to w'^^ = 
(equation pT]) ) is perfectly absorbing for this solution (it is satisfied identically by the 
outgoing solution), whereas the freezing- condition Bi corresponding to ^2^'' = is 
not, although the correction in ^2'*^'' falls off quite fast with boundary radius (~ r~^). 
More generally for arbitrary ^, Bl=i is perfectly absorbing whereas BL<i is not. 

3. 2. Initial data for numerical evolutions 

The initial data for our numerical evolutions are taken from the exact linearized 
multipolar wave solutions described in section 13.11 evaluated at t = 0. We shall 
consider the cases £ = 2, 3, 4. For definiteness, we choose m = 2 throughout; we have 
checked that our results are insensitive to the particular choice of m. 
We take an outgoing wave with mode function 



Font{x) = Aexp 



{x - ro) 



(60) 



The parameters are taken to be rg — 15 and a = 1.5. These choices are dictated by 
the requirement that at the outer boundary radius R of our computational domain, 
which we place at i? = 30, the wave initially vanish to numerical roundoff error; this 
guarantees that the initial data are consistent with the boundary data at i? = 30. 
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Figure 1. Outgoing linearized gravitational wave with multipolar indices £ = 
m = 2 and mode function l|60| l with parameters A = 1, ro = 15 and a = 1.5. 
Shown are the exact solutions I I57I I for the Zerilli scalar $(+) (solid) and I I58I I 
for the Regge- Wheeler scalar <E>(~) (dashed) as functions of time, evaluated at 



R = 30. 



Figure [T] shows the exact solutions for the RWZ scalars (for £ = 2) as functions of time 
evaluated at the outer boundary r = R = 30. 

Although we will evolve the full nonlinear Einstein equations, wc do not solve the 
full nonlinear constraints initially (of course the exact solution satisfies the constraints 
to linear order). This is because we are interested only in comparing the numerical 
solution with the exact solution to linear order in perturbation theory. 

3.3. Numerical evolution method 

Our numerical implementation uses the Caltech-CorncU Spectral Einstein Code 
(SpEC), which is based on a pseudo-spectral collocation method described in more 
detail for example in [24]. We evolve the first-order generalized harmonic form of the 
Einstein equations as described in [19] . The gauge-source functions H" in (|30p are 
taken to be zero when that equation is evaluated in Cartesian components, compatible 
with the TT gauge of the exact linearized solution. 

The computational domain for the evolutions presented here is a sphere of radius 
R = 30. It is subdivided into a small central sphere of radius Ar = 7.5 surrounded 
by three spherical shells, each of extent Ar. In the spherical shells, each Cartesian 
tensor component of each tensor field is expanded in Chebyshev polynomials in the 
radial direction and in spherical harmonics Yim{d, 4>) in the angular directions. In the 
central sphere, the numerical solution is expanded in the basis functions introduced by 
Matsushima and Marcus [33] with parameters a = 1 and (3 = 2: each Cartesian tensor 
component is expanded as f{r,9,(j)) = Y^^gj^ fnimQnd'^)Yim{0, 4>), where Qniir) 
depends on the index £ and Qniif) ^ r^ near the origin. The number of radial 
expansion coefficients is the same in each subdomain, and is denoted by N,.. Likewise, 
the highest retained spherical harmonic index is the same in each subdomain, and is 
denoted by iV^. In all subdomains we retain all spherical harmonic m coefficients with 
|m| < £. 

For the RWZ formalism, we need to compute the expansion coefficients of various 
tensors with respect to the basis Y^YatSa, Yab, Sab defined in (fT2|) . The SpEC code 
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is already capable of computing expansions in spin-weighted spherical harmonics. The 
conversion between the two bases is detailed in [Appendix A| We note that there is no 
need to use spin-weighted harmonics, it is simply convenient in the SpEC code. 

The evolution equations arc integrated in time using a fourth-order Runge-Kutta 
scheme, with a Courant factor Af/Axmin of at most 2.25, where Axmin is the smallest 
distance between two neighbouring collocation points. (This method of defining a 
Courant factor for a pscudospectral code with uneven grid spacing is only a crude 
approximation; so in practice, numerical stability can be achieved with Courant factors 
larger than unity.) As described in [21], the top four coefficients in the tensor spherical 
harmonic expansion of each of our evolved quantities are set to zero after each time 
step; this eliminates an instability associated with the inconsistent mixing of tensor 
spherical harmonics in our approach. The auxiliary ODEs ()10p are integrated using 
the same Runge-Kutta scheme. 

The BCs on the incoming fields of the generalized harmonic Einstein equations at 
the outer boundary and at the internal boundaries between neighbouring subdomains 
are enforced using a penalty method |34| . At the internal boundaries, the incoming 
fields on each boundary are set (weakly via the penalty method) equal to the 
corresponding outgoing field on the neighbouring boundary. At the outer boundary, we 
always impose the constraint-preserving boundary conditions (|38p on the constrained 
degrees of freedom, and we vary the BCs on the physical and gauge degrees of freedom 
as described below. 

3.4- Boundary conditions on gauge fields 

Besides the BCs on the physical degrees of freedom that are the subject of this paper, 
and those on the constraint degrees of freedom for which we choose equation ([55]) . 
we must also specify a BC PO)) on the gauge degrees of freedom. The multipolar 
gravitational waves described in section [3.11 are in TT gauge, 

Sgtt^O, <5.g-'=0, V^Sg'^^O, (61) 

where V is the flat-space covariant derivative. This implies that the harmonic 
gauge condition (|30p with vanishing gauge source functions H" — when evaluated in 
Cartesian components — is satisfied in linearized theory. However, the gauge boundary 
data F2 in (|40p do not vanish for this solution; we have to specify them explicitly. 
One might be tempted to replace the gauge BC (|40|) with 

i""a^ = 0, (62) 

which is compatible with TT gauge. Unfortunately, the ensuing IBVP is ill posed 
[35] . Alternatively, one might hope to transform the solution to a different gauge such 
that both the harmonic gauge condition and the homogeneous version of the gauge 
BC (|40]) are satisfied; however, it turns out that this cannot be achieved for a purely 
outgoing solution [55] . 

Hence we have decided to specify gauge boundary data that are computed 
from the exact linearized solution. Alternatively, we set = 0. The latter is 
compatible with the exact solution at early times provided the initial data have 
compact support and the outer boundary is sufficiently far out. Once the outgoing 
wave hits the outer boundary, there will be a (small) reflection which is pure gauge 
and hence should not affect the physical gravitational radiation as represented by the 
RWZ scalars. 
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3.5. Numerical results 

Here we describe the results of evolving initial data describing exact solutions of the 
linearized Einstein equations (section (|3.ip ) using our fully nonlinear evolution code 
that includes our new absorbing BCs. Because we will compare the RWZ scalars 
extracted from our nonlinear numerical evolution with their exact (linearized) values, 
our comparisons will disagree because of nonlinear effects. Therefore we repeat each 
run at several different amplitudes in order to separate these nonlinear effects from 
the boundary reflections of interest. The numerical resolution (Nr^Ni) is chosen 
high enough such that the numerical truncation error is negligible compared to the 
differences we are studying. 

We begin with the quadrupolar solution (£ ~ 2). Here we consider amplitudes 
in the range lO^'^ ""' ^ A ^ 10~^'^^. We performed a convergence test which showed 
that for any amplitude in this range, a numerical resolution of Nr = 81, A'^ = 8 is 
sufficient to remove the effects of numerical truncation error. 

First we consider B2, which was shown in section [3Tl to be perfectly absorbing 
for the £ = 2 solution. As mentioned in section 13.41 there is some freedom in the 
choice of gauge boundary data . First we compute these functions from the 
exact solution; with this choice of boundary gauge, the entire spacetime metric should 
agree with the exact linearized solution. We compute from the numerical evolution the 
£ = m = 2 component of the RWZ scalars, evaluated at the outer boundary, for both 
an odd-parity and an even-parity wave (this is the only nonvanishing component in the 
exact solution, up to the reality condition ()5ip ). We then compute the same quantity 
from the exact linearized solution, and plot the difference in figure [2l This difference 
has been normalized by the amplitude A, so that if the quantity plotted decreases 
at least linearly with decreasing amplitude, the difference between the numerical and 
exact solution decreases at least quadratically with amplitude, i.e. the two solutions 
agree to linear order in perturbation theory. This is indeed what we see in figure [2] 
(the decay exponent of the normalized difference is found to be close to 2). Similar 
behaviour is found for the remaining components of the RWZ scalars, which vanish 
for the exact solution. 

Next we replace the analytic gauge BC with the homogeneous one, F^ = 0, also 
referred to as a freezing gauge BC. The results are virtually unchanged (figure [3|). 
This demonstrates that the ambiguity in the choice of gauge BC has no effect on the 
gauge-invariant quantities, in all further evolutions discussed below, we will continue 
to use the freezing gauge BC. 

As discussed in section 13. 1[ we expect the freezing gauge BC to cause a small 
gauge wave reflection. This should show up in gauge-dependent quantities such as the 
spacetime metric. Figure |4] shows the component gtx of the metric, which vanishes 
for the exact solution because of the TT gauge. Indeed this quantity is found to 
decay quadratically with amplitude when analytical gauge BCs are used (left panel 
of figure m note again the quantity plotted has been normalized by the amplitude) . 
For the freezing gauge BC (right panel), the curves begin to overlap for sufficiently 
small amplitudes. This indicates that the numerical solution differs from the exact 
solution to leading (linear) order in perturbation theory. We had to use somewhat 
smaller amplitudes for this plot (10~^'^^ ^ ^ ^ 10"^ ''^) so that the nonlinear effects 
are smaller than the gauge wave reffection that we are looking for. 

We now repeat the evolution shown in figure [3l but using Bi; this corresponds to 
freezing at the boundary. As shown analytically in scction lSTSl Bi is not compatible 
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Figure 2. Difference between the numerically-extracted RWZ scalar <1?22 and 
its exact linearized value, evaluated on the outer boundary, and divided by the 
amplitude A. The exact linearized solution and the initial data for the evolution 
contain only the £ = 2, m = 2 mode of odd (left panel) or even (right panel) 
parity. The evolution uses B2 and analytic gauge BCs. Roundoff effects become 
visible at the smallest amplitude. 
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Figure 3. Same as figure|2] except freezing gauge BCs arc used in the evolution. 




Figure 4. The component gtx of the metric, normalized by the amplitude A, for 
the even-parity evolutions shown in figures [2] and |3] The two panels are identical 
except for the gauge BCs used: analytic (left panel) and freezing (right panel). 
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Figure 5. Same as figure |3] except the evolution uses BC Bi. 



with the exact outgoing solution. Figure [5] shows the RWZ scalars for this numerical 
evolution. The curves corresponding to smaller amplitudes nearly overlap, which 
implies that the numerical solution differs from the exact one at linear order. Note 
however that the normalized difference is still relatively small (about 10~^ or 10~^, 
depending on the parity). The remaining modes of the RWZ scalars not shown here 
behave in the same way as for B2, i.e. they decay at least quadratically with amplitude 
(note these modes vanish for the exact linearized solution). 

We now turn to evolutions of octupolar (£ ~ 3) waves. Here we consider 
amplitudes in the range 10~"^'^ ^ A ^ 10~^, and we find that a numerical resolution 
of Nr = 91, Ni = 10 is sufficient to remove the effects of numerical truncation error 
for the amplitudes in this range. Figure [6] shows evolutions using Si, B2, and B3. 
(We continue to use freezing gauge BCs.) As expected from the analytic discussion in 
section [3?T1 S3 is found to be consistent with the exact linearized solution, whereas Bi 
and B2 are not. To show the dependence on the order of the BC more clearly, in figure 
[7] we plot the time integrals of the curves in figure [H] versus the amplitude. We see 
clearly how for Bi and B2, this time integral saturates in the limit of small amplitudes, 
at a level which is lower for B2 than for Bi . This saturation to a horizontal line at small 
amplitudes again indicates a difference between the numerical and exact solution that 
decays only linearly with amplitude. The reason why these curves increase for larger 
amplitudes is because nonlinear effects become important. The curve for B3 behaves 
differently: it continues to decay with decreasing amplitude. (The slight levelling off of 
this curve at the smallest amplitude is caused by numerical round-off effects.) This is 
consistent with the theoretical prediction that B3 is perfectly absorbing for waves with 
^ ^ 3, so the difference between the numerical and exact linear solution is (almost) 
entirely due to nonlinear terms. Again, we find a nearly quadratic decay of the S3 
curve in figure [71 corresponding to a nearly cubic decay of the difference between the 
numerical and exact solution. 

Finally we evolve the £ = 4 solution. Here we consider amplitudes in the range 
10^"* ^ A ^ 10~^-^ and the numerical resolution is Nr = 101, N( ~ 12. Figure 
[8] compares Bi (corresponding to freezing ^Pq Sit the boundary) and S4 (which is 
predicted to be perfectly absorbing for this solution). Again the numerical results 
confirm that S4 is compatible with the exact linearized solution whereas Bi is not. 
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Figure 6. DifTerence between the numerically-extracted RWZ scalar <1>32 ^ and 
its exact linearized value, evaluated on the outer boundary, and divided by the 
amplitude A. The exact linearized solution and the initial data for the evolution 
contain only the £ = 3,m = 2 odd-parity mode. The different panels correspond 
to absorbing BCs Bj^ with L = 1, 2, and 3. Freezing gauge BCs arc used. 




Figure 7. Time integral of the curves shown in figure \E\ plotted against the 
amplitude A. The different curves correspond to absorbing BCs Bl with L = 1, 
2 and 3. 
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Figure 8. Same as figure (6] except the exact linearized solution and the initial 
data for the evolution now contain only the odd-parity £ = A,m = 2 mode. The 
different panels correspond to Bj^ with L = 1 and 4. 



3.6. Comparison with the predicted reflection coefficients 

In [1], solutions to the linearized Bianchi equations with higher-order absorbing BCs 
equivalent to Bl were studied theoretically. For monochromatic radiation, analytical 
expressions for the reflection coefflcients relating the amplitude of an outgoing wave to 
the spurious reflected ingoing wave were derived (equation 96 of [1]). In this section, 
we compare these predicted reflection coefficients with the numerical results presented 
in section [531 

The predicted reflection coefficients in [T] are given as functions of 27ri?/A, where 
A is the wavelength of the radiation and R is the radius of a spherical outer boundary. 
Our numerical tests, however, evolve wave packets and not monochromatic waves. So 
to obtain the reflection coefficients as a function of A for our numerical evolutions, 
we proceed as follows. We measure the radiation reflected off the outer boundary by 
taking the difference between our numerically-evolved solution and the exact linearized 
solution at the boundary, which we have plotted as a function of time in flgurcs [51 
[51 and [HI In order to minimize the nonlinear effects, the lowest amplitudes are used 
{A = 10-2-89 for ^ = 2, A = 10-3-^ for £ = 3 and A = IQ-^ for i = 4). Taking 
the Fourier transform of any of the curves in these figures (prior to taking absolute 
values) yields the difference as a function of frequency, from which we can obtain 
the difference as a function of wavelength, A<I>(A). Similarly, we can obtain $(A) by 
taking the Fourier transform of the exact linearized solution ^(t), evaluated at the 
boundary. Given these quantities, we define the reflection coefficient at any particular 
wavelength A as A<i>(A)/$(A). 

The values of both the analytical and numerical reffection coefficients are plotted 
as a function of X/R in figure [9l for radiation with angular momentum numbers 
£ = 2,3,4 and absorbing BCs BL<e- This demonstrates the importance of higher- 
order BCs when the numerical simulation contains multipolar radiation. For example, 
assume the outer boundary is located at twice the characteristic wavelength of the 
simulation, so that X/R. = 0.5. Then, for radiation with £ = 3 and freezing- BCs 
{Bi), the predicted errors due to spurious reflections off the boundary are 3 x 10"^. 
The use of B2 decreases these errors 100-fold, down to 3 x 10"^. For multipolar 
radiation with f = 4, a factor of about 3600 in accuracy is gained as one goes from Bi 
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Figure 9. Reflection cocflicients as a function of \/R (where _R is the outer 
boundary radius and A is the wavelength) for odd-parity linearized solutions to the 
Einstein equations, with absorbing BCs Bl of selected orders L. Curves denote 
coefficients qi^ i calculated from equation 96 of [T], and points denote coefficients 
computed from numerical evolutions of lincarized-wavc initial data. Different 
panels correspond to solutions with different values of the angular momentum 
number £. Note that B^—i is perfectly absorbing and hence Ql i = 0. 

(where the error is 9 x 10"'*) to B3 (where the error is 2.5 x 10~^). 

Figure [9] also shows that the predicted and numerical reflection coefficients agree 
well except for small values oi X/R and small values of the reflection coefficient. The 
discrepancies are likely to be caused by a combination of different effects, (i) We 
have computed the numerical reflection coefficients from a fully nonlinear evolution 
whereas the predictions are only valid in linearized theory. These nonlinear effects 
should become less important as the amplitude of the wave is decreased, (ii) Numerical 
roundoff error becomes important for small values of the difference between the exact 
and numerical solution (visible in figures [51 [H and [S]) . (iii) There is an accumulated 
error of spurious reflections passing through the origin and reflecting again, which 
is not accounted for in the theoretical reflection coefficients. One can see this in 
the secondary pulse at t ~ 45 in figures [5l [6] and [51 (iv) For Af/(27ri?) > 1, the 
theoretical predictions for the refiection coefficients start to break down because the 
ingoing solutions used for the calculation in [T] are not valid as the wave approaches 
r = 0. 
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4. Conclusions 

An algorithm for numerically implementing the higher-order absorbing BCs Bl for 
Einstein's equations presented in [1] has been defined, tested, and shown to work. 
Our method is based on a reformulation of the BCs in terms of the gauge-invariant 
RWZ scalars. This approach relies on the assumption that close to the outer boundary, 
the Einstein equations can be linearized about a given background spacetime, which in 
this paper is taken to be flat. It also requires a spherical outer boundary. A key feature 
of our algorithm is the introduction of auxiliary variables intrinsic to the boundary, 
a technique familiar from the fields of, for example, computational electromagnetism 
and acoustics. 

We have used a generalized harmonic formulation of the Einstein equations |19j : 
however, it should be straightforward to adapt our algorithm to different formulations, 
with minor modifications. Similarly, our use of spectral numerical methods is not an 
essential ingredient. Our boundary algorithm can be added to any existing numerical 
relativity code in a modular way: all that needs to be done is to implement the 
auxiliary evolution system at the boundary and provide for the necessary exchange 
of information with the main evolution system. In order to estimate the additional 
computational cost in doing so, we note that the number of auxiliary variables 
that need to be evolved for the absorbing EC Bl (which is perfectly absorbing for all 
multipoles ^ ^ i) is 2L{L + 1)^. (The ranges are 1 s$ fc ^ L, < £ ^ L, ^ m ^ £.) 
E.g. for i = 3, this amounts to 96 functions (of time only), which is small compared 
to a typical number of grid points in the interior domain. 

We have implemented our method in the Caltech-Cornell SpEC code and have 
tested it by evolving initial data with not only quadrupolar {i = 2), but also higher 
multipolar {£ = 3,4) radiation, imposing the hierarchy of BCs Bl with L < i. We 
use linearized solutions of the Einstein equations [211 HI] as initial data, and evolve 
these in a fully nonlinear code. We demonstrate that with decreasing amplitude of the 
radiation in the initial data, and perfectly absorbing BCs BL=i, the numerical solution 
decays at least quadratically to the exact linearized solution, as expected. For Bl<i^ 
the difference between the numerical and exact solutions decays only linearly as the 
amplitude is decreased, indicating that linear, spurious reflections are being introduced 
into the solution. We have estimated the magnitude of these spurious reflections in 
our numerical simulations and find good agreement with the theoretical reflection 
coefficients derived in [1], within the range of validity of the two calculations. For 
multipolar radiation with ^ = 3 or 4, we have seen that even without imposing a 
perfectly absorbing BC Bl=i, one can decrease the reflections off the outer boundary 
dramatically by increasing the order L of the BC Bl- For instance, for an outer 
boundary radius of twice the wavelength, R = 2A, and ^ = 3, we achieve a 100-fold 
decrease in reflection by using B2 rather than Bi. 

The most important application of our work is to the simulation of isolated 
systems, such as coalescing binary black holes. In these simulations, the most 
sophisticated BC on the gravitational radiation currently in use is the freezing-\['o 
BC [23l[24l[25l[26l[271[28l[29], which corresponds to Bi. Note that the ratio of 
the wavelength to the outer boundary radius in the numerical tests presented in 
this article is X/R ~ 0.1; such a small value for X/R makes it easy to ensure that 
the BCs are consistent with the initial data because the radiation vanishes at the 
outer boundary. In contrast, for binary black hole simulations, X/R is much closer to 
unity, in particular during the inspiral phase. One can see from the graphs in figure 
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ini that errors due to reflections with imperfectly absorbing BCs (including Bi) are 
several orders of magnitude larger at, say, X/R ~ 0.5 than they are at X/R ~ 0.1. 
Furthermore, the majority of current binary black hole simulations use far cruder 
BCs on the gravitational radiation than Bi: Typically, a sequence of adaptive mesh 
refinement [37] grids of decreasing resolution are used that extend out to larger and 
larger radii (see e.g. [38l [39l |40l IHl |42] ) . On the coarsest grid, the waves are no longer 
properly resolved, and some ad hoc BC (such as a Sommcrfeld BC) is imposed. It is 
clear that this approach will cause spurious reflections, and given the results of [27], 
we suspect that these will be considerably larger than those caused by Bi. 

Inspiralling and merging equal mass non-spinning binary black holes emit 
predominantly quadrupolar {£ = 2) radiation. For a typical outer boundary location 
used in numerical simulations of these events, R = 2A, the errors due to spurious 
reflections of quadrupolar radiation with Bi are predicted in [1] to be 6 x 10^^, which 
seems very small. However, these reflections may interact with the evolution in such 
a way as to result in an accumulating phase error. For instance, the runs in |28| 
showed phase errors of a few hundredths of a radian which were attributed to the 
outer boundary location. Moreover, as the mass ratio M1/M2 deviates from unity, 
the energy emitted in non-quadrupolar modes rapidly increases: while for M1/M2 — 1, 
less than 0.1 per cent of the energy is emitted with £ > 2, this fraction increases to a 
few per cent for M1/M2 = 2 and exceeds 10 per cent for the still fairly modest mass 
ratio M1/M2 = 4 [43l|44]. Because Bi has larger reflection coefficients at a given 
value of A/i? for £ > 2 modes than for £ = 2 modes (cf. figure [9]), this BC will not 
perform as well for unequal mass binary black hole simulations; in these cases, the 
higher-order BCs presented in this paper will be even more important than they are 
for equal mass simulations. Finally, when considering the merit of the higher-order 
BCs on simulations containing different modes (£, m), one needs to remember that the 
wavelength of a mode depends on the value of m. During the inspiral, modes with 
smaller \m\ have larger wavelengths A (typically in proportion to 1/|to|), and therefore 
the ratio of wavelength to boundary radius, X/R, depends on m. From figure [9] we 
then deduce that modes with the same £ but different m will have different reflection 
coefficients. For example, the (2, 1) mode has been shown to contribute significantly to 
the unequal mass binary black hole inspiral radiation [43]. Because its wavelength is 
longer than for the (2, 2) mode, the reflection coefficient for Bi will be correspondingly 
higher. 

For these reasons, we are confident that higher-order BCs will improve the 
accuracy of long-term unequal mass binary black hole simulations. Whether such 
an improvement is desirable will of course depend on the precision requirements for 
the computed waveforms. Parameter estimation for LIGO requires phase errors of 
hundredths of a radian [45], and LISA will require yet more accurate waveforms. 

Our work can also be applied to the problem of Cauchy-perturbative matching. 
In this approach, one matches solutions to the full nonlinear Einstein equations on 
an interior spatial domain to solutions of the linearized equations on an exterior 
domain extending out to large radii. These linearized equations can be represented for 
instance by the gauge- invariant RWZ equations. The question now arises as to how 
one should impose outer BCs on the Einstein equations in the interior with boundary 
data computed from the exterior linearized solution. Sections 12.41 and 12.51 contain a 
specific prescription for how this can be done. 

In future work, we plan to generalize our algorithm by (i) including first-order 
corrections in M/R for the curvature and backscatter on a Schwarzschild background 
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spacetimc (of mass Af), (ii) allowing for arbitrary coordinates of the background 
spacetime, and (iii) adapting the algorithm to more general (i.e., not necessarily 
spherical) boundary shapes. The generalized RWZ formalism |18| upon which our 
calculations arc based and the absorbing BCs formulated in [J [5] allow for such 
generalizations. In more realistic situations where an exact solution is not at hand, 
one can assess the quality of the BCs by comparing with a numerically computed 
(fully nonlinear) reference solution on a large domain |27| . Ultimately, wc plan to 
apply our implementation of absorbing BCs to binary black hole simulations using 
the Caltech-Cornell SpEC code. 
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Appendix A. Conversion to spin-weighted harmonics 

The spin- weighted spherical harmonics |46j sYf,m{0,4') can be constructed from the 



standard spherical harmonics recursively using the relations 

oYlm = Ylra, (A.l) 

= \{l-s){l^s^\)\-^^''^sYi^, (A.2) 

=-[(^+s)(^-S+l)]-l/2e,y,„. (A.3) 

The operators (3 and 9 arc defined by 

SsYf,™ = {-de - i CSC 6*90 ^ scold) sYim, (A. 4) 

BsYfm = {-de + i CSC 6*90 - s cot 6')sY«m. (A. 5) 
Next we set up a tetrad e,, = (t, r, m, m). (Greek indices /i, i^, . . . from the middle 



of the alphabet denote tetrad indices.) In components with respect to spherical polar 
coordinates (t, r, 0, 0), 

i„ = (1,0,0,0), r„ = (0,1,0,0), m„ = -^(0,0,1, isinfl), (A.6) 

and m is the complex conjugate of m. An arbitrary rank-2 tensor T can be expanded 
as 

T= ^ TZ:" si.,,.^)Yira^^.®^.. (A.7) 
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Here the spin weight is defined by s(e^,ei,) = s(e^) + .s(e^) where s(m) = — 1, 
s(m) = +1, and s(t) = s(r) = 0. The SpEC code provides routines that compute the 
expansion coefiicients T^,^^" in (jA.7p for a given tensor T. 
The basis harmonics defined in (|12p can be written as 



YA,hn = 7= [lYernmA - -lYimmA), (A. 8) 

V2r 

SA/m ^ -^-^^^-^^^^{iyemrhA + -lYimmA), (A.9) 

V2r 

YABJm = ^^^2r'2^ ^ ^-{2YlmmA'mB + ^2Yl>mmAmB), (A. 10) 



SAB,hn = 7;-^ {2YhnmAmB - -2y«mTOAms). (A. 11) 



2r2 

(Recall that A = (-^ — 1)(£ + 2).) We also have gAB = 2r~^m(^77iB). Using these 
relations it is straightforward to express the amplitudes of the perturbation (fT3 l) -([T4 l) 
and p9 |) -([2T |) in terms of the expansion coefficients &g1^" , 

/i*,,™ = ^=^^(^5*™ + ^^tn.)^ (A.12) 

K,,rr. = + <5.g[™) , (A. 13) 



ir 


ir 




V2£(^4 


-1) 


— ir^ 




2VA£(^4 


-1) 



i?tt,£m = (A. 15) 

Htr,lm=59ta, (A.16) 

i?rr,^m = ^ffl™, (A.17) 

Qt,,^ = -j^=={59\f, ~ <55^*™), (A. 18) 

Q.,,™ = -^^^(J^rm _ (^.19) 

= ^.Cf , (A.20) 

G,,„ = {5gtr + ■ (A.21) 
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